Cite this Article
Li Chen-Hui, Han Xiu-Jun, Luan Ying-Wei, Li Jian-Guo. Abnormal breakdown of Stokes–Einstein relation in liquid aluminium. Chinese Physics B, 2017, 26(1): 016102  
Permissions
Abnormal breakdown of Stokes–Einstein relation in liquid aluminium
Li Chen-Hui, Han Xiu-Jun †, Luan Ying-Wei, Li Jian-Guo
Laboratory of Advanced Materials Solidification, School of Materials Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

 

† Corresponding author. E-mail: xjhan@sjtu.edu.cn

Abstract

We present the results of systematic molecular dynamics simulations of pure aluminium melt with a well-accepted embedded atom potential. The structure and dynamics were calculated over a wide temperature range, and the calculated results (including the pair correlation function, self-diffusion coefficient, and viscosity) agree well with the available experimental observations. The calculated data were used to examine the Stokes–Einstein relation (SER). The results indicate that the SER begins to break down at a temperature (∼1090 K) which is well above the equilibrium melting point (912.5 K). This high-temperature breakdown is confirmed by the evolution of dynamics heterogeneity, which is characterised by the non-Gaussian parameter α . The maximum value of α , α , increases at an accelerating rate as the temperature falls below . The development of α was found to be related to the liquid structure change evidenced by local five-fold symmetry. Accordingly, we suggest that this high-temperature breakdown of SER has a structural origin. The results of this study are expected to make researchers reconsider the applicability of SER and promote greater understanding of the relationship between dynamics and structure.

PACS: 66.20.Cy;66.10.cg;02.70.Ns;61.20.Ja
Keyword:diffusion;viscosity;dynamics heterogeneity;local five-fold symmetry
1. Introduction

Dynamics in liquid metals, such as diffusion and viscosity, play a critical role in governing crystal nucleation and growth, as well as glass transition.[1, 2] In simple liquids, the self-diffusion coefficient D is associated with the viscosity η through the Stokes–Einstein relation (SER) ,[3, 4] where is the Boltzmann constant; d is the effective particle diameter, which is assumed to be a temperature-independent constant; and c is a constant whose value depends on the boundary conditions: 2 for slip boundary conditions and 3 for stick boundary conditions. SER was first proposed to describe the Brownian motion of suspended particles that are spherical in form and have a diameter of d.[5, 6] The relation has been widely tested in application to many different types of systems, such as molecular liquids,[7] Lennard–Jones liquids,[8] hard-sphere fluids,[9] and metallic glass-forming melts.[10] In many cases, SER applies quite well at temperatures above the critical temperature ( ,[11, 12] where is the glass transition temperature) predicted by the mode coupling theory (MCT).[13] Because is always far lower than the melting point , it is generally believed that SER describes the relationship between D and η well at . However, previously, we found that SER breaks down for Cu Zr melt at the critical temperature , which is well above ,[14] and such an abnormal breakdown has also been observed in another multicomponent glass-forming metallic liquid.[15] The abnormal breakdown reported in Refs. [14] and [15] has been ascribed to enhanced dynamics heterogeneity. Because enhanced dynamics heterogeneity is present in a wide variety of liquid systems, we wondered whether a similar breakdown could be detected in a simpler system, such as a monatomic melt. A survey of the literature suggests that the validity of SER for monatomic systems has rarely been tested. We undertook this study because we believed that an exploration of the applicability of SER to a monatomic metallic melt such as that which occurs in liquid aluminium (Al) would have a universal and fundamental physical significance.

As a popular engineering material, pure Al has been widely studied as a prototype metal.[16–22] A precise determination of D and η is required to test the validity of SER for liquid Al. Recently, the D and η of liquid Al have been measured by several researchers.[23–26] It should be stressed that these researchers’ measurements were limited to narrow temperature ranges and did not involve the undercooled state. We are not aware of any experimental data on D and η for undercooled liquid Al in the literature. There are two main aspects to the difficulty of obtaining experimental data on undercooled liquid metal. First, the metastable undercooled state is hard to maintain for long periods of time during experiments. Second, considerable uncertainty exists in the experimental data because of convection inside the sample in terrestrial measurements. Molecular dynamics (MD) simulation was employed in this study to consider a broad temperature range and avoid experimental difficulties. The advantage of MD is the determination of the trajectory of all atoms in a simulation system, which is difficult to achieve by other methods. MD simulations have provided in-depth knowledge of the structural and transport properties of liquid metals, and recently the MD method has garnered considerable attention.

Although ab initio MD simulation can achieve high accuracy, the system size (usually ∼200 atoms) that can be simulated at present is too small to solve the dynamics and structure of liquid metal. The results obtained will inevitably exhibit a size effect. Therefore, classical MD simulation with a well-accepted Finnis–Sinclair potential was performed in this study, using a large simulation cell containing more than 30,000 atoms.

2. Simulation methods

All of the simulations were carried out using the large-scale atomistic/molecular massively parallel simulator (LAMMPS) code.[27] Details of the potential used in the simulations are provided in [28]. Based on our previous work, the of Al for this potential was taken to be 912.5 K.[29] The simulation box (8.1×8.1×8.1 nm) had 32000 atoms, with periodic boundary conditions in three directions. The isothermal–isobaric (NPT) ensemble was used. The external pressure was maintained at zero using the Parrinello–Rahman algorithm.[30] The temperature was controlled by the Nosé–Hoover thermostat.[31, 32] A time step of 1 fs was applied, using the velocity Verlet algorithm[33] to integrate Newton’s motion equation. The system was heated to 2050 K and held at 2050 K for 0.6 ns to obtain a homogeneous liquid. The liquid was then cooled to 50 K at a cooling rate of K/s. During the cooling, atomic configurations were recorded at intervals of 20 K for use in further dynamic and structural analyses. The recorded configurations were then equilibrated at the corresponding temperature for 0.3 ns to obtain the equilibrium state. During the equilibration and sampling processes, the volume and energy were monitored and found to remain constant, i.e., the aging effect was negligible (see Appendix A). The structure-like pair correlation function was derived and averaged from ten equilibrium configurations. The dynamics-like mean square displacement was computed by averaging over a large number of trajectory snapshots. We used three independent samples of a system of 32000 atoms to improve accuracy and minimize errors. The evaluation time was carefully modified, based on the temperature considered, to obtain a time-independent value. The final calculated values of D and η were the averages for the three independent samples.

3. Results and discussion
3.1. Narrative description of the cooling

Figure 1 shows the profile of the system during cooling. The profile is characterised by the energy per atom ( ) and the density (ρ) of the system. As the temperature decreases, E decreases and ρ increases. Discontinuous changes in E and ρ were observed at 600 K, which means that the system cannot maintain a liquid state at that temperature, and crystallisation takes place. Because of the pure circumstance in MD simulation, substantial undercooling of 300 K ( ) is normal. It should be noted that the cooling rate ( K/s) that we employed is close to the upper limit that can be achieved by MD simulation. Because of the limitations of computer power and computing time, a cooling rate between and K/s is widely used. However, configurations recorded during the cooling process at a slower cooling rate will require a shorter time to reach the equilibrium state, which is important in the subsequent calculation of relevant ensemble-averaged quantities.

Fig. 1 (color online) Total energy per atom ( ) and density (ρ) are plotted as a function of temperature.
3.2. Structure
3.2.1. Pair correlation function

The pair correlation function g(r) has been used to describe the structural feature of the liquid system. It is expressed as follows:

(1)
where g(r) is the probability of finding an atom in the range from r to N is the number of atoms; ρ is the average density of the system; R i and R j are the positions of atoms i and j, respectively; and δ is the Dirac delta function. For a monatomic system, g(r) can be simplified as follows:
(2)
where is the number of atoms within the distance between r and . As a contrast, we show the simulated and experimental g(r) values at 943 K, 1023 K, and 1323 K in Fig. 2. Overall good agreement is observed between the simulated and experimental g(r) values,[34, 35] which indicates that our simulation approach is reasonable.

Fig. 2 (color online) The simulated g(r) and the experimental one for liquid Al. The data are sparsely depicted to clearly show the comparison. The experimental data for 943 K and 1323 K are from [34], for 1023 K from [35].

Figure 3 illustrates the variation in g(r) and how the peaks and shoulders develop from high temperature to low temperature. As Figure 3 shows, the height of the first peak in g(r) increases step-by-step with decreasing temperature, implying an enhancement of short-range order in the first coordination shell. A gradually enhanced shoulder is found on the right-hand side of the second peak, indicating a gradual change in the liquid structure as the liquid passes from superheating to supercooling. This gradual change is generally considered to be associated with the formation of atomic clusters.[36, 37] We will return to this point later.

Fig. 3 (color online) Pair correlation functions g(r) at diverse temperatures (as shown by the arrow, the temperature changes from 1950 K to 650 K with a step of 100 K). The inset is the calculated coordination number.

The inset figure illustrates the coordination number (CN) as a function of the temperature. CN is derived from g(r) using the following equation:

(3)
where n is the atomic number density and is the first minimum position in g(r). As the inset figure shows, CN increases linearly from 10.68 ± 0.10 to 12.69 ± 0.10 with decreasing temperature as the liquid changes from superheated at 1950 K to supercooled at 650 K.

3.2.2. Angle distribution analysis

The change in the liquid structure with cooling involves not only the number and length of the bonds but also their directions. To obtain a comprehensive picture of the structural change, the dihedral-angle distribution and the bond-angle distribution were studied. As Figure 4(a) shows, the dihedral-angle distribution function has a distinct peak at approximately 60° and a second peak at approximately 108°. As the temperature decreases, the first peak shifts to a larger angle and the second peak shifts to a smaller angle, indicating a continuous change in the liquid structure. As Figure 4(b) shows, the bond-angle distribution function has two peaks, centred at approximately 56.30° and approximately 107.60°. It should be noted that the peak positions in Fig. 4(b) are not far from the characteristic bond angles of icosahedrons (namely, 63.5° and , suggesting the possible presence of a large number of icosahedral-like clusters. As the figure shows, those peaks become higher with decreasing temperature, which demonstrates the strengthening of clustering. In addition, a gradually enhanced hump is observed at approximately 150°. The hump is pronounced for the undercooled liquid but spreads out for the superheated liquid. Therefore, this hump is believed to be related to the structural character of the undercooled liquid.

Fig. 4 (color online) Evolution of the angle distribution functions for liquid Al with decreasing temperature. (a) The dihedral-angle distribution, and (b) the bond-angle distribution.
3.3. Dynamics
3.3.1. Mean square displacement and self-diffusion coefficient

The time dependence of the mean square displacement (MSD) reflects important information on atomic diffusivity. The self-diffusion coefficient D can be calculated from the long-time evolution of MSD[38] as follows:

(4)
using the Einstein equation:
(5)
where is the ensemble average.

For a solid system, MSD saturates to a finite small value, while if the system is liquid, MSD increases significantly with time. As shown in Fig. 5(a), the significant increase in MSD indicates that the systems, including the deeply undercooled ones, were in a liquid state during the calculation of the relevant ensemble-averaged quantities. Because of the high cooling rate and the pure circumstance during simulation, liquid Al can remain in the deep undercooled regime for a long time. For short time intervals, MSD increases in proportion to t 2, which is typical of the ballistic motion of atoms. For long time intervals, MSD exhibits a linear dependence on t, i.e., the atom motion becomes diffusive. The calculated values of D are shown in Fig. 5(b). The experimental results[23, 24] are shown for comparison. The calculated results are in good agreement with the experimental results. Both the calculated and experimental data obey the Arrhenius equation. Because of the broad temperature range considered in the simulation, the Arrhenius fitting of the calculated D values can be separated into two parts intersecting near a temperature of 1090 K.

Fig. 5 (color online) (a) Log-log plot of the MSDs for liquid Al at varying temperatures. (b) Self-diffusion coefficients ( ) extracted from MSD plotted in the logarithmic scale against the inverse temperature. The squares are the present simulated results. The triangles and circles are the experimental data of Refs. [23] and [24], respectively. The lines are their Arrhenius fits.
3.3.2. Viscosity

The viscosity η is calculated by integrating the autocorrelation function of the pressure tensor, based on the Green–Kubo relation[39]

(6)
where is the off-diagonal element ( , α, , y, z) of the pressure tensor at time t, V is the volume of the system, and is the Boltzmann constant. Using Eq. (6), three values of η can be determined independently from the three off-diagonal pressure components (P xy , P xz , and . In real liquids, η is considered to be isotropic, so these three independent η must converge to one single value. A dimensionless quantity ξ is introduced to estimate the convergence error quantitatively:
(7)
where
(8)
η xy is the viscosity calculated from P xy , η xz is the viscosity calculated from P xz , and η yz is the viscosity calculated from P yz . As an example, we present the three normalised autocorrelation functions of the pressure tensor at 650 K in Fig. 6(a). Because 650 K is the lowest temperature considered in this study, it required the longest integration time for the autocorrelation function to decay. As Figure 6(a) shows, the integration time was sufficiently long to characterise the decay process, and the differences among the three pressure components were negligible. The calculated η and ξ for 650 K are plotted against the integration time in Fig. 6(b). The figure shows that η eventually reaches a constant value. The convergence error ξ is 4.75%. Figure 7 shows the average normalised autocorrelation functions (i.e., the arithmetic means of , , and ) at different temperatures. Figure 7 shows that as temperature increases, the decay of the autocorrelation functions accelerates.

Fig. 6 (color online) (a) The three normalised autocorrelation functions of pressure tensor ( , , and ) at 650 K. (b) Viscosities η calculated from the three off-diagonal pressure components, and the convergence error ξ are plotted as a function of the integration time for liquid Al at 650 K.
Fig. 7 (color online) The average normalised autocorrelation functions of pressure tensor at diverse temperatures.

Figure 8(a) shows the convergence errors at different temperatures. The convergence error ξ fluctuates within a small range and decreases with increasing temperature. The calculated viscosities are presented in Fig. 8(b), along with the experimental data.[25, 26] The calculated viscosities agree well with the experimental data and fall between the two sets of experimental data. The calculated viscosities can be fitted to the MCT power law:

(9)
This fitting yields parameter values of and K. The temperature dependence of η can be described by fitting the η values to the Vogel–Fulcher–Tammann (VFT) law[40–42]
(10)
This fitting yields parameter values of B = 1314 K and K. Both the MCT fitting and the VFT fitting describe the temperature dependence of the calculated viscosities well over the entire temperature range considered in the simulation. Following the usual convention, we fitted η utilising the Arrhenius equation. Because the Arrhenius fitting is highly dependent on the fitting range, we tried many times to modify the fitting. We ultimately found that one Arrhenius expression was not sufficient to characterise the evolution of η. Generally speaking, the Arrhenius fitting fails to capture the evolution of η in the low-temperature range. In the deeply undercooled state ( K), the evolution of η is clearly different from that indicated by the high-temperature Arrhenius fitting (HTAF). To be more specific, upon cooling at K, η increases at a much more rapid rate than that described by the HTAF. The change from Arrhenius to non-Arrhenius dependence indicates a crossover from uncorrelated to correlated dynamics in the liquid. Because of the variation in the calculated values, it is difficult to determine the crossover temperature directly by examining the discrepancies between the calculated η values and the HTAF-predicted values. However, we noted that the MCT and VFT fittings yield very reasonable predictions in the low-temperature range. We considered it to be better to use the degree of deviation of the MCT or VFT curve from the HTAF line to determine the crossover temperature. However, because the HTAF line intersects the MCT curve, it is not convenient to use the MCT curve for this purpose. We therefore defined the onset deviation point as the temperature at which the discrepancy between the HTAF line and the VFT curve was up to 2%. Based on this definition, the onset deviation temperature (i.e., the crossover temperature was found to be ∼1090 K.

Fig. 8 (color online) (a) The convergence error (ξ) versus temperature. (b) The computed viscosities using the Green–Kubo relation are plotted against temperature, and their temperature dependence is fitted by the MCT, VFT, and Arrhenius expressions. The error bar is to show the discrepancy in the calculated value among three independent samples, and the same below.
3.3.3. Self-intermediate scattering function and α-relaxation time

Figure 9(a) shows the temperature dependence of the self-intermediate scattering function , which is expressed as

(11)
characterises the mean relaxation time of a system. The wave number q in is the position of the first peak in the static structure factor S(q). As shown in Fig. 9(a), the time needed for to decay increases with decreasing temperature because of weakened atomic mobility. Over short periods of time, exhibits a Gaussian-like decay representing the ballistic motions of atoms. Over long periods of time, decays to zero in nearly exponential fashion on a picosecond timescale, suggesting that the motion of atoms becomes diffusive. The slow diffusion stage is also called α-relaxation. The α-relaxation time τ can be empirically defined as the time required for to reduce to .[43, 44] The increase in τ with decreasing temperature is illustrated in Fig. 9(b). As with η, calculated values of τ can be fitted using the MCT, VFT, and Arrhenius expressions. We do not attempt to present a detailed discussion of the fitting of τ here because it is similar to that of η. However, we emphasise that with decreasing temperature, the deviation of τ from the HTAF line becomes very clear, and the failure of the Arrhenius equation is demonstrated again. Using the same criterion as in the analysis of η described above, we found that for τ, the crossover from Arrhenius to non-Arrhenius behaviour also occurs at ∼1090 K. It should also be noted that, according to MCT, τ is related to according to the following equation:
(12)
Regression analysis yielded parameter values of and K for this equation. The fitted γ value is nearly the same as that derived from the calculated viscosities. The fitted is close to but slightly higher than the value extracted from η. It is worth noting that the fitted value is lower than the crystallisation temperature of 600 K. According to the fitted value, SER should hold quite well within the temperature range considered in this study, as we discuss in the next section.

Fig. 9 (color online) (a) Self-intermediate scattering function at different temperatures, and the time is plotted in the logarithmic scale. The α-relaxation time τ is determined by the cut of , as illustrated by the dashed line. (b) α-relaxation time τ versus inverse temperature, and the lines are the MCT, VFT, and Arrhenius fittings.
3.4. Breakdown of Stokes–Einstein relation

Once D and η have been calculated, the validity of SER can be evaluated. According to SER, under slip boundary conditions, the effective diameter d can be determined from the equation . As long as d is a temperature-independent constant, SER is thought to hold; otherwise, SER breaks down. As Figure 10(a) shows, the calculated d is close to the first nearest-neighbour distance in g(r). At high temperatures, d fluctuates around a value of approximately 2.85 Å. At low temperatures, d deviates significantly from the constant. More precisely, d exhibits a decrease at (∼1090 K). To assess the possibility of a size effect in the MD calculation, we conducted additional tests for systems of 4,000 and 108,000 atoms. There was no obvious size effect evident in the results for d (see Appendix B). The reduction in d is indicative of increasingly correlated atom motions. Because η is proportional to τ via the Maxwell relation,[2, 15] it is reasonable to assess the validity of SER using the parameter , which is plotted with respect to temperature in Fig. 10(b). Clearly, the evolution of with temperature can be divided into two regimes with different characteristic slopes. At (∼1090 K), increases with decreasing temperature at a constant rate that is four times larger than that at . Based on the results of d and , the onset of SER breakdown is considered to occur at . It must be emphasized that the crossover temperature is higher than and by approximately 550 and 177 K, respectively. To the best of our knowledge, this is the first time that the abnormal breakdown of SER in a monatomic liquid metal has been observed.

Fig. 10 (color online) (a) Effective Stokes–Einstein diameter d determined by the calculated D and η. The dashed line is the first nearest-neighbour distance derived from g(r), 2.8 Å, which is nearly independent of temperature. The dotted-dashed line is the fitting result for d. Correlation coefficient for , and 0.79 for . (b) Parameter versus temperature. The lines are the linear fittings with and 0.86 for the low and high temperature ranges separated by , respectively.
3.5. Dynamics heterogeneity and clustering
3.5.1. Non-Gaussian parameter

In a previous study, we found that the abnormal breakdown of SER occurred in Cu Zr melt as a result of a sudden increase in dynamics heterogeneity.[14] We subsequently investigated dynamics heterogeneity in liquid Al using the non-Gaussian parameter :[45]

(13)
where is the MSD and is the mean quartic displacement. For an isotropic system, . When , a number of atoms have a higher mobility than others, or the atom mobility varies. The variation in mobility is illustrated by a colour coding of atom displacement over a time interval of 1000 fs at 650 K in Fig. 11(a). Figure 11(a) displays a broad distribution of atom displacements, suggesting varying atom mobilities.

Fig. 11 (color online) (a) Snapshot of atom position at 650 K. Colours represent the displacement of atom during the time interval of 1000 fs. (b) Non-Gaussian parameters α for liquid Al versus time at varying temperatures. (c) Maximum of the non-Gaussian parameter α depicted in the semi-logarithmic scale against . The red line is the Arrhenius fitting for high temperature (HT fit), and blue line for low temperature (LT fit). At temperature around , a kink is observed.

Like those of MSD and , the evolution of can be divided into different regimes based on its time dependence. Initially, starts at zero because of the ballistic motion of atoms and increases slightly because of the anisotropies of ballistic motion. then increases with time according to a power-law relation (proportionally to ). Eventually, decays to be approximately proportional to and tends to zero because of diffusive motions.[46] Therefore, between the short-time ballistic regime and the long-time diffusive regime, there exists a maximum value of , referred to here as . In general, increases with decreasing temperature, and the time required to reach becomes longer. provides information on correlated jumps of atoms over time and complicated clusters in space.

As Figure 11(b) shows, the evolution of is consistent with the above descriptions. As the temperature decreases, increases by an order of magnitude, which is evidence that the dynamics of liquid Al becomes progressively more heterogeneous. To illustrate its temperature dependence, is plotted with respect to temperature in Fig. 11(c). As the figure shows, increases upon cooling at a considerably enhanced rate at (∼1090 K), and SER begins to fail. Atoms with high mobility make large contributions to D, whereas η and τ are strongly influenced by atoms with low mobility. Therefore, an enhanced dynamic heterogeneity can result in the breakdown of SER.

3.5.2. Icosahedral cluster and local five-fold symmetry

It is widely accepted that liquid metals are composed of atomic clusters.[47–52] Atoms in a cluster usually have lower mobility because they have to coordinate with each other. Consequently, dynamics heterogeneity in liquid metal is believed to be initiated by the formation of atomic clusters. Icosahedral clusters, which are an important and fundamental type of atomic cluster, are thought to be responsible for the undercooling ability of liquid metal.[53] It is therefore important to study icosahedral clusters.[54–58] Research has shown that icosahedral clusters play an important role in the dynamics–structure relationship.[59–63] For example, Jakse et al. noted that an increase in icosahedral clusters is responsible for non-Arrhenius dynamics slowing down in liquid and undercooled Cu55Hf45 and Cu62Hf38 alloys.[62] Similarly, Cheng et al. demonstrated that the icosahedral cluster type could be considered the structural origin of non-Arrhenius dynamics.[63] Because non-Arrhenius dynamics is closely related to the breakdown of SER, it is important to examine the possible correlation between the breakdown of SER and icosahedral clusters. Therefore, the focus of the analysis described below was the icosahedral cluster type.

As the structural results for g(r) and the angle distribution suggest, atomic clusters in liquid Al gradually become pronounced as the temperature decreases. To obtain more detailed information, we characterised atomic clusters by means of a Voronoi tessellation analysis.[64, 65] A Voronoi polyhedron is similar to a crystallographic Wigner–Seitz cell. The Voronoi index is a type of coordinate polyhedron surrounding a central atom, where n i is the number of the face with i vertices. The total number of faces is equal to the CN of the central atom. The Voronoi index of a perfect icosahedral cluster is . Voronoi polyhedra with indices of , , , , and have been identified as icosahedral-like clusters (or distorted icosahedra).[66–68] As Figure 12(a) shows, perfect icosahedra ( ) only occupy a small proportion of liquid Al. The majority of these Voronoi polyhedra are distorted, and distorted icosahedra marked by are the most frequent. We can see from Fig. 12(a) that , , , and increase with decreasing temperature. Unlike these four types, and are nearly temperature independent.

Fig. 12 (color online) (a) Fractions of icosahedral and icosahedral-like clusters at different temperatures. (b) Evolution of LFFS with decreasing temperature.

To reflect various clusters in a unified way, a structural indicator of local five-fold symmetry (LFFS)[69] can be defined as follows on the basis of Voronoi tessellation analysis:

(14)
According to Eq. (14), LFFS is the ratio of the number of pentagons to the number of Voronoi faces. LFFS has been proven to be sensitive to glass transition and mechanical deformation.[70–73] Figure 12(b) illustrates the evolution of LFFS with cooling. The figure shows that LFFS increases with decreasing temperature and that the temperature dependence of LFFS has two linear regimes: a high-temperature regime and a low-temperature regime. Within the high-temperature regime, the slope is , whereas in the low-temperature regime, the slope is . The two regimes intersect at a temperature of ∼1090 K, which is the same temperature as . Based on the previously described relationship between atomic clusters and dynamic heterogeneity, it is reasonable to conclude that a correlation exists between the development of LFFS and dynamic heterogeneity. This correlation suggests that the abnormal breakdown of SER has a deeper origin which is linked to change in the structure of the liquid.

4. Conclusion

We studied the structure and dynamics of liquid Al using molecular dynamics simulation. The calculated results for the structure included the pair correlation function, the dihedral-angle distribution, and the bond-angle distribution. The structural results suggest that the liquid structure changes significantly from high temperature to low temperature. The dynamics results included the self-diffusion coefficient D, the viscosity η, and the α-relaxation time τ. The calculated D and η values were in good agreement with experimentally measured values. We used the calculated D, η, and τ values to assess the validity of the Stokes–Einstein relation, SER. The results reveal that SER breaks down at a critical temperature K, which is higher than as determined by MCT fitting and is also higher than the melting-point temperature . Such a high temperature breakdown of SER has rarely been reported for a monatomic liquid metal.

To explain this abnormal breakdown, we used the non-Gaussian parameter to describe the dynamics heterogeneity. The results show that dynamics heterogeneity is accelerated with decreasing temperature at . Accordingly, this abnormal breakdown of SER is considered to be rooted in enhanced dynamics heterogeneity. Because dynamics heterogeneity is closely associated with atomic clustering, we investigated icosahedral and icosahedral-like clusters by means of Voronoi tessellation analysis. We found that the number of clusters increases considerably with cooling. LFFS, which is defined as the ratio of the number of pentagons to the number of Voronoi faces, was introduced to facilitate a uniform description. The results suggest that the evolution of LFFS is consistent with the development of dynamics heterogeneity. Based on this finding, we propose that change in the structure of the liquid is the intrinsic reason for the abnormal breakdown of SER.

Reference
[1] Glicksman M E 2011 Principles of Solidification New York Springer
[2] Ngai K L 2011 Relaxation and Diffusion in Complex Systems New York Springer
[3] Tyrrell H J V Harris K R 1984 Diffusion in Liquids London Butterworths
[4] Balucani U Zoppi M 1994 Dynamics of the Liquid State Oxford Clarendon
[5] Einstein A 1956 Investigations on the Theory of Brownian Movement New York Dover
[6] Sutherland W 1905 Philos. Mag. 9 781
[7] Chang I Sillescu H 1997 J. Phys. Chem. B 101 8794
[8] Bordat P Affouard F Descamps M Müller-Plathe F 2003 J. Phys. -Condens. Matter 15 5397
[9] Kumar S K Szamel G Douglas J F 2006 J. Chem. Phys. 124 214501
[10] Chathoth S M Samwer K 2010 Appl. Phys. Lett. 97 221910
[11] Cicerone M T Ediger M D 1996 J. Chem. Phys. 104 7210
[12] Gotze W Sjogren L 1992 Rep. Prog. Phys. 55 241
[13] Das S P 2004 Rev. Mod. Phys. 76 785
[14] Han X J Schober H R 2011 Phys. Rev. B 83 224201
[15] Jaiswal A Egami T Zhang Y 2015 Phys. Rev. B 91 134204
[16] Li F Liu X J Hou H Y Chen G Chen G L 2011 J. Appl. Phys. 110 013519
[17] Shen B Liu C Y Jia Y Yue G Q Ke F S Zhao H B Chen L Y Wang S Y Wang C Z Ho K M 2014 J. Non-Cryst. Solids 383 13
[18] Jakse N Pasturel A 2013 Sci. Rep. 3 3135
[19] Recoules V Crocombette J P 2005 Phys. Rev. B 72 104202
[20] Li H Wang G H Zhao J Bian X F 2002 J. Chem. Phys. 116 24
[21] Smith P M Elmer J W Gallegos G F 1999 Scripta Mater. 40 937
[22] Cherne F J III Deymier P A 2001 Scripta Mater. 45 985
[23] Demmel F Szubrin D Pilgrim W C Morkel C 2011 Phys. Rev. B 84 014307
[24] Kargl F Weis H Unruh T Meyer A 2012 J. Phys.-Conf. Ser. 340 012077
[25] Wang D Overfelt R A 2002 Int. J. Thermophys. 23 1063
[26] Yamasaki T Kanatani S Ogino Y Inoue A 1993 J. Non-Cryst. Solids 156-158 441
[27] Plimpton S 1995 J. Comput. Phys. 117 1
[28] Mendelev M I Kramer M J Becker C A Asta M 2008 Philos. Mag. 88 1723
[29] Li C H Han X J Luan Y W Li J G 2015 Chin. Phys. B 24 116101
[30] Parrinello M Rahman A 1980 Phys. Rev. Lett. 45 1196
[31] Nosé S 1984 J. Chem. Phys. 81 511
[32] Hoover W G 1985 Phys. Rev. A 31 1695
[33] Swope W C Andersen H C Berens P H Wilson K 1982 J. Chem. Phys. 76 637
[34] Waseda Y 1980 The Structure of Non-Crystalline Materials New York McGraw-Hill
[35] Morris J R Dahlborg U Calvo-Dahlborg M 2007 J. Non-Cryst. Solids 353 3444
[36] Schenk T Holland-Moritz D Simonet V Bellissent R Herlach D M 2002 Phys. Rev. Lett. 89 075507
[37] Chen K Y Liu H B Li X P Han Q Y Hu Z Q 1995 J. Phys. -Condens. Matter 7 2379
[38] Rapaport D C 2004 The Art of Molecular Dynamics Simulation 2004 Cambridge Cambridge University Press
[39] Allen M P Tildesley D J 1987 Computer Simulation of Liquids Oxford Clarendon
[40] Vogel H 1921 Phys. Z. 22 645
[41] Fulcher G S 1925 J. Am. Ceram. Soc. 8 339
[42] Tammann G Hesse G 1926 Z. Anorg. Allg. Chem. 156 245
[43] Egelstaff P 1967 An Introduction to the Liquid State London Academic Press
[44] Hansen J P McDonald I R 2006 Theory of Simple Liquids 3 London Academic Press
[45] Rahman A 1964 Phys. Rev. A 136 405
[46] Caprion D Matsui J Schober H R 2000 Phys. Rev. Lett. 85 4293
[47] Lee G W Gangopadhyay A K Kelton K F Hyers R W Rathz T J Rogers J R Robinson D S 2004 Phys. Rev. Lett. 93 037802
[48] Kelton K F Lee G W Gangopadhyay A K Hyers R W Rathz T J Rogers J R Robinson M B Robinson D S 2003 Phys. Rev. Lett. 90 195504
[49] Liu H R Liu R S Zhang A L Hou Z Y Wang X Tian Z A 2007 Chin. Phys. 16 3747
[50] Zhang J X Li H Zhang J Song X G Bian X F 2009 Chin. Phys. B 18 4949
[51] Hou Z Y Liu R S Tian Z A Wang J G 2011 Chin. Phys. B 20 066102
[52] Wang L Zhang Y N Mao X M Peng C X 2007 Chin. Phys. Lett. 24 2319
[53] Frank F C 1952 Proc. R. Soc. A 215 43
[54] Kim T H Kelton K F 2007 J. Chem. Phys. 126 054513
[55] Ganesh P Widom M 2006 Phys. Rev. B 74 134205
[56] Leocmach M Tanaka H 2012 Nat. Commun. 3 974
[57] Shen Y T Kim T H Gangopadhyay A K Kelton K F 2009 Phys. Rev. Lett. 102 057801
[58] Di Cicco A Trapananti A Faggioni S 2003 Phys. Rev. Lett. 91 135505
[59] Han X J Li J G Schober H R 2016 J. Chem. Phys. 144 124505
[60] Cheng Y Q Ma E 2011 Prog. Mater. Sci. 56 379
[61] Pasturel A Tasci E S Sluiter M H F Jaske N 2010 Phys. Rev. B 81 140202
[62] Jakse N Nguyen T L T Pasturel A 2013 J. Appl. Phys. 14 063514
[63] Cheng Y Q Sheng H W Ma E 2008 Phys. Rev. B 78 014207
[64] Finney J L 1970 Proc. R. Soc. A 319 495
[65] Finney J L 1977 Nature 266 309
[66] Delogu F 2009 Phys. Rev. B 79 064205
[67] Itoh K Hashi K Aoki K Mori K Sugiyama M Fukunaga T 2007 J. Alloys Compd. 434-435 180
[68] Hao S G Wang C Z Kramer M J Ho K M 2010 J. Appl. Phys. 107 053511
[69] Peng H L Li M Z Wang W H 2013 Appl. Phys. Lett. 102 131908
[70] Cheng Y Q Ma E Sheng H W 2008 Appl. Phys. Lett. 93 111913
[71] Peng H L Li M Z Wang W H 2011 Phys. Rev. Lett. 106 135503
[72] Wakeda M Shibutani Y Ogata S Park J 2007 Intermetallics 15 139
[73] Cheng Y Q Cao A J Sheng H W Ma E 2008 Acta Mater. 56 5263